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Abstract 

A subspace method for channel estimation has been recently proposed III for tackling the pilot 
contamination effect, which is regarded by some researchers as a bottleneck in massive MIMO systems. 
It was shown in III that if the power ratio between the desired signal and interference is kept above 
a certain value, the received signal spectrum splits into signal and interference eigenvalues, namely, 
the “pilot contamination” effect can be completely eliminated. However, III assumes an independently 
distributed (i.d.) channel, which is actually not much the case in practice. Considering this, a more 
sensible finite-dimensional physical channel model (i.e., a finite scattering environment, where signals 
impinge on the base station (BS) from a finite number of angles of arrival (AoA)) is employed in this 
paper. Via asymptotic spectral analysis, it is demonstrated that, compared with the i.d. channel, the 
physical channel imposes a penalty in the form of an increased power ratio between the useful signal 
and the interference. Furthermore, we demonstrate an interesting “antenna saturation” effect, i.e., when 
the number of the BS antennas approaches infinity, the performance under the physical channel with 
P AoAs is limited by and nearly the same as the performance under the i.d. channel with P receive 
antennas. 
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I. Introduction 

Massive MIMO systems that employ a large number of antennas at the base station (BS) 
have attraeted signifieant interest reeently [|2l, flUl, dH. The main advantage of using massive 
MIMO lies in the signifieant improvement of speetral and energy effieieney. However, when we 
shift our attention to the multieell seenario, the massive MIMO system would, unfortunately, 
be plagued by the so-ealled “pilot eontamination” effeet, whieh is due to the use of non- 
fully orthogonal pilot sequenees aeross all the eells. However, it is often hard to aehieve full 
orthogonality of the pilot sequenees among the terminals aeross the eells, due to the limited 
eoherenee time of the mobile eommunieation ehannels. The performanee of massive MIMO 
systems under the full reuse of pilot sequenees was studied intensively by Marzetta [|2l, where 
it is shown that, by linear proeessing of the reeeived signal at the BS, the performanee is 
only interferenee-limited and does not depend on the transmitted power of users, i.e., signal-to- 
interferenee ratio (SIR) eannot grow unboundedly. 

Several motivating works have been eondueted aiming to address the above pilot eontamina- 
tion problem. In [|6l, authors use a eoordinated ehannel estimation seheme whieh exploits the 
information embedded in the seeond-order statisties (e.g., eovarianee matrix) of uplink ehannels. 
The main idea of the seheme in |[6l is to assign pilots to users assoeiated with eovarianee matriees 
exhibiting maximum orthogonality of signal subspaees. Although this seheme greatly alleviates 
the pilot eontamination, but it ineurs too mueh eoordination sinee all eovarianee matriees of all 
uplink ehannels must be learned beforehand. 

In 0, by leveraging a key observation that the quasi orthogonality among the ehannel veetors 
implies that the ehannel veetors of the desired users are eigenveetors of the eovarianee matrix 
of the reeeived signal in the asymptotie limit, Ngo and Larsson 0 proposed a blind ehannel 
estimation method whieh mitigates the need of pilots. However, it is worth noting this seheme 
heavily hinges on a large number of BS antennas as well as large bloek length (i.e., sample 
size), therefore its performanee might be unsatisfaetory in the “not very large” regime. 

Muller et. al. [|8l proposed another blind pilot deeontamination method whieh, in essenee, aims 
to distinguish users in the amplitude domain by exploiting the differenee between the ehannel 
gains of the intended users and the interfering users. Therefore, this work ean be regarded as 
a parallel to |l6l, whieh essentially attempts to distinguish users with the same pilot sequenee 
in the angular domain by exploiting the non-isotropy of angles of arrival (AoA) multipaths of 
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the channel. By means of random matrix theory (RMT) and free probability theory (FPT), (Sl 
demonstrates that under a mild power ratio the subspace-based estimation scheme is capable of 
totally removing the interference and hence completely tackling the pilot contamination problem. 

We note that the channel considered in [l8]| consists of independently distributed (i.d.) Rayleigh- 
faded entries, implicitly implying a rich scattering environment. However, in the real world, it 
is quite common that the number of scatterers is limited and correspondingly the number of 
AoAs is finite [|9l|, Q, ifT^ . Therefore, considering a more sensible physical channel model (or 
alternatively, called finite-dimensional channel, due to the finiteness of its degrees of freedom) 
is of significant importance since the finiteness of degrees of freedom of the channel is relevant 
to the ability of subspace method to identify the eigenvalues of desired and interfering users. In 
the absence of a priori knowledge of the channel at the BS, it is natural to assume the AoAs in 
the channel model to be uniformly distributed in the interval [0, tt]. Therefore, the total number 
of AoAs, which roughly corresponds to the total number of scatterers around the BS, will be 
the only key parameter (besides the ratio of the channel gains that correspond to the desired 
users and the interfering users, respectively) that might have crucial impact on the performance 
of the blind subspace method. Hence a natural question one may ask: What is the impact of 
finite-dimensional channel on the performance of the subspace-based channel estimation scheme? 
What is the price paid for this finite dimensionality of the channel? Moreover, does increasing 
the number of antennas at the BS help to alleviate the performance degradation? 

To answer the above questions, the core task is to characterize the spectrum of the observed 
signal under the physical channel model. Despite its difficulty, we manage that by leveraging 
tools in RMT and FPT. Intuitively, the case of the i.d. channel can be regarded as a special case 
of our result, i.e., when the number of AoAs approaches infinity. The contributions of this paper 
are summarized as follows: 

1) It is shown that to guarantee the performance of the subspace-based channel estimation 
scheme under the physical channel, we should pay a cost of an extra power margin between 
the intended users and the interfering users with respect to (w.r.t.) a given BS. 

2) For multicell multiuser MIMO system where the users of each cell are seen from distinct 
AoAs w.r.t. a given BS, it is shown that the performance is mainly determined by the cell 
associated with smaller number of AoAs. 

3) It is demonstrated that there exists an antenna saturation effect at the BS, i.e. adding antennas 
beyond a threshold at the BS is of limited help to enhance the performance under the physical 
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channel model. To be eonerete, the performanee under a physieal ehannel with P AoAs is 
nearly the same as the performanee under the i.d ehannel model with P antennas at the BS 
(see FiglU and SecUV] for a detailed description of the simulation setup and discussion). 



Fig. 1. Antenna saturation effect under physical channel. L = 4 cells, K = 5 users per cell, P = 200 AoAs, coherence time 
N = 400 symbol periods and per-user SNR = —5dB. ZF is used for channel estimation and MF for data detection. Array 
elements are critically-spaced. 



Fig. 2. The effect of the number of AoAs on the required power ratio Pj/Ps- M — 200 receive antennas, L — 2 cells, K = 5 
users per cell, coherence time N = 400 symbol periods, per-user SNR = OdB. ZF is used for channel estimation and MF for 
data detection. Array elements are critically-spaced. 


In the following we provide some intuitive remarks of the above results. When the number of 
AoAs is deereased, the eorrelation (or more exaetly, the eoherenee) among the ehannel veetors 
will inerease eorrespondingly. As a result, the condition number (or eigenvalue spread) of the 
ehannel matrix Hg (eorresponding to inside-eell users) as well as Hj (eorresponding to outside- 
eell users) will both inerease. Therefore, the left and right endpoints of the eigenvalue elusters 
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corresponding to Hg and Hj, respectively, tend to be closer, even overlapping (e.g., see FigjS]). 
In order to keep the above two eigenvalue clusters apart, we have to pay a cost of boosting 
the channel gain ratio between the desired users and the interfering users. Figl2l illustrates the 
BER performance of the subspace method proposed in |l8l when we vary the number of AoAs, 
denoted P. It can be seen from Figl2]that the smaller P is, the higher channel gain ratio Pg/Pi 
is required. 

As a final remark, it is worthwhile of pointing out the different role that channel correlation 
plays in the two pilot decontamination methods, namely, the Bayesian channel estimation in Q 
and the subspace method in |[8l as well as this paper, that is, while the correlation is beneficial in 
the former method, it is however, disadvantageous in the latter. The key reason for this difference 
lies in that, in a nutshell, the former is a Bayesian method while the latter is a non-Bayesian one 
(i.e., no a priori assumed). Specifically, the correlation is advantageous for the linear (MMSE) 
channel estimation (since it helps to distinguish users in the angular domain) at the cost of 
acquiring the a priori (namely, the covariance matrices) and coordinating the pilots. On the other 
hand, as a non-Bayesian method in essence, the subspace method performs the channel estimation 
in a nonlinear way, which highly relies on the instantaneous property (such as the eigenvalue 
spread of the channel matrix), rather than the statistical property in the Bayesian method 
Since correlation tends to increase the eigenvalue spread, it will impact the ability of subspace 
method to distinguish the desired users and interfering users via eigenvalues clustering. 

The rest of the paper is organized as follows: In Sec HI] we introduce the physical channel 
model and review the subspace-based channel estimation. In Secjllll the asymptotic eigenvalue 
distribution (AED) of the channel is derived in the Stieltjes domain. This will enable us to derive 
analytical expressions, from which the support of the distribution can be identified, hence yielding 
the main results of this paper. Eurther, SecjTVj leverages the obtained formulae to demonstrate 
the impact of the physical channel on the spectral spread of signal and interference subspaces, 
and BER simulation results are presented for performance comparison. SecjVj summarizes the 
main results and concludes this paper. 

II. PROBLEM FORMULATION 

A. Channel model 

We will consider the uplink in multicell multiuser MIMO communication system with L cells. 
Each BS is equipped with a uniform linear array with a large number of antennas M, serving K 
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single-antenna users. We also assume time-division duplex arehiteeture. The ehannel is assumed 
narrowband flat-fading, whieh remains eonstant over a eoherenee time of N symbol periods, 
and changes independently from one coherence time to another (i.e., block fading channel ifTTI l. 
Furthermore, we assume time-synchronous users, and full reuse of pilot sequences among all 
cells to facilitate channel estimation at the BS. 

Since the channel is estimated blindly, we assume no a priori information about the channel 
is available at BSs, including the angular spread, which is in contrast to the Bayesian channel 
estimation method O. Thus, we assume all AoAs are uniformly distributed in [0,7r], which is 
a reasonable, yet mathematically tractable assumption. However, the theoretical analysis in our 
paper gives insights into the performance of the subspace method under a physical channel with 
small angular spreads as well. 

The array response in a given direction is quantified by the so-called steering vector. Hence, 
the channel from the kth user in the ith cell to the intended BS, denoted h^i, can be modeled 
as a linear combination of all steering vectors (see e.g., ifTOll . Il6l ) 



( 1 ) 


where Pki is the number of i.i.d. AoA multipaths, s{ipkij) G is the steering vector correspond¬ 
ing to the angle of arrival ipkij associated with the jth path, 1/y/Pki serves as a normalization 
factor and akij ~ CAA(0,/3fcj) denotes the channel gain associated with the jth direction, where 
\/%^i is the average channel attenuation, due to path loss and shadowing effect. Moreover, it 
is assumed that all Lfkij and akij are independent over user index k, cell index i and direction 
index j. The length-M steering vector associated with the angle of arrival i.pkij is given by 

• • •, where we use “j” to denote the imaginary 

unit and fi{^pkij) is a function of Lfkij. 


Let 



( 2 ) 


be the M x K effective fast-fading channel from K users in the ith cell to the intended BS, 
where Ski = {s^Lpkn), • • •, s^ifkip^j))/\/^i G comprises all steering vectors of the kth 

user in the ith cell w.r.t. the intended BS, and hki G consists of i.i.d. CAA(0,1) entries. 
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Fig. 3. A simple illustration of multipath environment where the signals of users 1 and 2 arrive at the BS from identical AoAs 
due to joint scatterers. 


The M X N signal received by the Zth BS during N consecutive symbol intervals may be written 
as 

L 

Yi = J2 (3) 

i=\ 

where Z?* = diag(/3ii,/92i, • • •,G comprises large-scale fading coefficients, Ti = 

diag(pij, • • • ^PKi) G consists of the average transmission powers, and Xi G is the 

normalized input symbols (preceded by short pilot signals), all for K users in the Zth cell. The 
entries of the noise matrix W) are assumed i.i.d. CAA(0,1). 

For the sake of analytical simplicity, we shall assume that all users belonging to the same (say, 
Zth) cell are seen from the same set of directions (whose cardinality is Pi, and a corresponding 
steering matrix Si) w.r.t. the intended BS. Therefore, in this distinct AoAs scenario, ^ takes the 
form of Hi = SiHi, where Hi G consists of i.i.d. CAA(0,1) entries, correspondingly, ([3]) 

can now be rewritten as: 

L 

Yi = J2 SAD^^Ti^/^X, + Wi. (4) 

i=l 

As a first step in addressing the channel model (SI), weTl make a further assumption that all 
Si are identical w.r.t. to the intended BS, which means all users in the network are seen from 
the same set of P AoAs by the intended BS (for the sake of illustration, see Fig©. Hence, in 
this identical AoAs scenario, (SI) can take a simpler form as: 

L 

Yi = sY, H^D]'^Ti^/^Xi + Wi. (5) 

i=l 

In SecHnl we will first focus on the channel model ([5]) and for the general case, i.e., the channel 
model (SI), we will treat it afterwards. 
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B. Subspace-based channel estimation 

In massive MIMO, by using subspace-based estimation approach, each user’s channel can be 
estimated blindly and accurately up to a scalar ambiguity, specifically, by identifying the signal 
subspace from interference subspace. This is possible because when the number of receive 
antennas is very large and the uncorrelated channel is assumed, the signals of different users are 
projected onto quasi-orthogonal subspaces. Another key condition that makes it possible to split 
the spectrum into signal and interference eigenvalue clusters, is due to the fact that inside-cell 
users and outside-cell users exhibit difference in the received power at the intended BS (e.g., 
due to path loss). These observations have been recently employed in massive MIMO [|3, |[8l 
to circumvent the problem of pilot contamination. 

Applying singular value decomposition (SVD) to the received signal © (see |[8l for more 
details), the first K eigenvectors, denoted Us G can then be identified, which correspond 

to the signal subspace. These eigenvectors represent an estimate of the channel up to scaling 
ambiguity which can be resolved by exploiting a short training sequence (e.g., 1 pilot) sent by 
each user. Having identified Us, the received signal is then projected onto Us, by which most 
of interference is annihilated and most of the thermal noise is removed. In this process, a new 
channel model is obtained: 

Yi = UlYi = GiXi + Wi ( 6 ) 

where Gi G and Wi G are the subspace channel and subspace noise matrix, 

respectively. Note that the original channel Hi is not needed for the detection of Xi, rather it 
can be estimated after Xi being detected for the purpose of downlink precoding, for example. 

III. Performance analysis 

In this section we will be primarily concerned with deriving the asymptotic eigenvalue dis¬ 
tribution (AED) of the channel models @ and ® with the aid of RMT and FPT. In fact, the 
density of the distribution is not explicitly derived. Instead, a fixed-point equation that satisfies 
the Stieltjes transfoml] ifT^ . [[T3ll of the density is given, which will serve as grounds for 
the derivation of the formula that helps identify the support of the distribution, in particular, the 


'Also known as the Cauchy transform, which encodes all the moments of the underlying distribution in a polynomial function. 
For a distribution function f{x), its Stieltjes transform is defined by G(s) = ^ 
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gap between the desired signal and interference subspaces. In Sec JTTT-AI we shall start with the 
channel model dS]) and the channel model dll) is treated in Sec JTTT-Rl 

A. Physical channel model with identical AoAs 

First, we should stress that the steering matrix S' in (O) is random, but fixed, which is dependent 
on the physical environment only. Equation dS]) can be written in the compact form as 

Yi = SHD^/^X + Wi (7) 

where H = (^i, .ffa, • • •, e X = (Xi, X 2 , • • •, e and D e 

j^KLxKL jg ^ diagonal matrix, where the components of the row vectors {pul3ii, • • •, PKil3Ki)A = 
1, 2, • • •, L are placed on its diagonal. 

If the true data covariance matrix S is assumed available at the BS, then the eigenvalue 
distribution of the channel can be accurately determined, for instance by S = {SH)D{SHy + 
Im, but in practice, one only has access to the sample covariance matrix computed from 
N data samples, which serves as an approximation of S. Thus, our objective is to study the 
eigenvalue distribution of Sa? in the asymptotic sense, where Sat is given by the following 
advanced information-plus-noise model: 

Stv = {SHD^^^X + Wi){SHb^/‘^X + Wi)^/N (8) 

which is difficult to handle by using Stieltjes transform approach ifTdll . ifTHl . 

Remark 1. The analysis of the eigenvalue distribution ofT^j^ under an arbitrary diagonal matrix 
D does not admit a tractable solution. Note that our primary concern here is to investigate the 
impact of power difference between the desired signal and the interference which guarantees, 
along with other system parameters, the separability of the signal and the interference subspaces. 
Therefore, without loss of generality, we assume the worst case scenario of power imbalance 
between inside-cell and outside-cell users. Specifically, we set all interference powers to the 
maximum interference power, denoted Pi, among all interfering users. On the other hand, we 
set the powers of all desired signals to the minimum signal power, denoted Ps, among all inside¬ 
cell users. Based on the above assumption, the diagonal matrix D is now comprised of two 
distinct masses, namely Pg and Pj, with multiplicity of K and K{L — 1), respectively. 

1) Assumption 
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For the sake of analytieal traetability and asymptotic results, we shall make the following 
assumptions: 

- The input matrix X consists of i.i.d. Gaussian entries. This assumption is not strict since 
the asymptotic result still holds given the entries are independent with finite moments of 
order greater than 2 ifT^ Section 5.1], which is fulfilled for most signal constellations. 

- The diagonal matrix D has bounded eigenvalues, and in the large limit, the eigenvalue dis¬ 
tribution converges to a deterministic distribution. Note that this is true since the attenuation 
due to path loss and shadowing is generally constant over many symbol intervals. 

- The entries of matrix S are assumed independent. Note that this assumption is violated in 
practice, where it is shown in [fTSll that the moments of Vandermonde matrix are bounded 
above and below by the moments of Marcenko-Pasture distribution [fT^ and Poisson dis¬ 
tribution, respectively. Nevertheless, when the number of receive antennas grows large, 
this assumption becomes reasonable flU, and hence the moments of Marcenko-Pasture 
distribution will become a good approximation of the moments of S. 

2) One-sided spectral analysis 

While the exact solution of the AED of ([8]) is difficult to obtain, to simplify the derivation 
we assume high SNR regime, i.e., W] = 0. Further, we also assume that K P and P is 
sufficiently large so that orthogonality between users’ channel vectors can still holcj^ Under 
these assumptions, the eigenvalue distributions of the signal and the interference can be studied 
separately. 

Proposition 1. Let S G Hi G and Xi G be defined as in (|7]). Also let 

M, K, P, N ^ oo with K P, K/M —)■ a, P/M ^ fi, and K/N 7. Consider the 
N X N Hermitian matrix F=Ps{SHiXiY{SHiXi)/MN. Then the AED of F converges to a 
non-random distribution with Stieltjes transform Gf{s) satisfying the following equation: 

(dy^iySG f{^s) -f 1) -f PsGf{s)[sGf{s) — y ~\~ l)(crsGjr(s) 

-fa — y){asG f{s) + a — fiy) = 0. (9) 


^It is shown in El that when M —> oo, spatial correlation yields only a minor penalty on the orthogonality condition compared 
to the independently distributed channel. Note also for two channels hki and hij as defined in (TJ with identical steering matrix 
S, Ymip^M^aah\^hij/M —>■ 0 almost surely sls K P. 
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Proof: consider the sealed matrix product 

Dk = y/ak/mpnAkBkCk ( 10 ) 

and let 

Gi = DlDk = akCl{AkBk)\AkBk)Ck/mpn ( 11 ) 

where is a multiplieative eonstant, and the matrices Ak G Bk G and Ck G 

are mutually independent. Moreover, the entries of Ak are assumed independent with zero mean 
and unit varianee, whereas the entries of Bk and Ck are assumed i.i.d. CAA(0,1). We also assume 
rank (AkBkCk) = h, i-C-, Gi has Ik nonzero eigenvalues. 

Following note that tr{Cl{AkBky{AkBk)Ck) = tr{{AkBky{AkBk)CkCl), i.e., both 
produet faetors inside the trace operators share the same nonzero eigenvalues and differ only in 
\n — lk\ zeros. Thus it would be more adequate to study the eigenvalue distribution of the Ik x Ik 
matrix 

G 2 = OkiAkBk^ iAkBk){CkCl)/mpn ( 12 ) 

without the need to obtain the eigenvalue distribution of Gi directly. Now let G 2 = G 21 G 22 , 
where G 21 and G 22 are defined, respectively, as 

G 21 = OkiAkBkY (AkBk) /mp, (13) 

G 22 = CkCl/n. (14) 

To find the AED of G 2 , one may first derive the marginal distributions of G 21 and G 22 by 
means of RMT, then FPT, namely, the free multiplicative law ifTTll . lIT^ . [fTSll ean be invoked 
straightforwardly. The usefulness of FPT lies in the possibility of computing the AED of a 
produet and a sum of Hermitian random matrices based solely on the marginal distributions, 
given the matriees are asymptotieally free (see e.g., [[I3l See. 2.4] for details on asymptotie 
freeness of random matriees). In FPT, the free multiplicative law of a product of random matrices 
is eomputed via the so-ealled S'-transform (see e.g., [fTSll . Def. 2.15.). More specifieally, when 
two Hermitian random matrices are asymptotieally free, with each matrix has an eigenvalue 
distribution that eonverges almost surely in distribution to a eompactly supported distribution, 
then the S'-transform of their product is the produet of their eorresponding S'-transforms [[US, 
Thm. 4.7]. 
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Under the assumptions made on and Bk, the asymptotie freeness condition is fulfilled (see 
Appendix for details). From ||9l, the Stieltjes transform Gg 2 i(^) of tho AED of G 21 satisfies the 
following equation: 

+ afcQ;is(Q;2 + 1 — 2ai)GQ^^{s) (15) 

+ {a2S + ak{ai - l)(ai - a2))GG2i(s) - 02 = 0 

where ai = Ik/m, and 0:2 = p/m are limiting ratios. Further, from Def. (21)j^ if we 
substitute s = —1/A in (fTSl) . and replace Gg2i(~ 1/-^) by “'^(^('^) + 1)’ we obtain 

afccr^A(T(A) + 1)^ + afcQ;iA(T(A) + l)^(cr2 — 2cri + 1) 

— A ( 0 : 2 /A — afc(Q'i — a2)(cri — 1))(T(A) + 1) + q;2 = 0. (16) 


Replacing A by T ^{z), and using BH Def. (20)], it is easy to check that the S'-transform Sg 2 i{^) 
reads 


ak{aiz + a2){aiz + 1)' 

Next, since the entries of C^/y/n are assumed i.i.d. each with zero mean and variance 1/n, 
then G 22 is a central Wishart matrix with n degrees of freedom BUI, BT^ . BTSll . Further, its 
eigenvalue distribution follows the well-known Marcenko-Pastur law BT^ . Thus when n,lk ^ 00 
with a fixed ratio Ik/n a^, the S'-transform Sg 22 {^) corresponding to G 22 reads as BEl, ifTSll 


Sg22{z) = 


(18) 


( 0 : 3 ^ -(-1) 

Returning to (fT^ . it is easy to see that, as a consequence of asymptotic freeness of G 21 
and G 22 (the validity of this assumption is discussed in details in Appendix), the S'-transform 
corresponding to G 2 is the product of S'g 2 i(^) and S'g 22 (^)> that is BT2l Thm. 4.7] 


<S'g2(^) ak{aiz + a2){a3Z + l){aiz + 1)' ^ 

and from BT^ Thm. 2.32], it can be shown by straightforward algebra the S'-transform corre¬ 
sponding to Gi (fTH) is 


SgAz) 


(^2 

+ Oi2){^z + 1)(2; -f as) ‘ 


( 20 ) 


^Note that different definition of Stieltjes transform is used in (9). 
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Now we are left with the task of determining the Stieltjes transform eorresponding to 
From dll Def. (20) and (21)], it ean be easily verified that the S'-transform Sp{z) and the Stieltjes 
transform Gf{s) of a distribution function F satisfies 

Gf{s) 


Sf {-sGf{s) - 1) = 


( 21 ) 


sGf{s) + 1 

Using (l20l) and with the aid of (|2TI) . it follows straightforwardly after some mathematical 
manipulations that the Stieltjes transform Ggi(s) corresponding to Gi satisfies the following 
equation: 

0:203(1 + sGgi) + ttkGciii — 03 + ’SGGj(ai 


O3 + OisGgi)(oi — O2O3 + OisGgi) — 0 


( 22 ) 


Finally, we remark that if we replace ak, oi, 02 and 03 in (l22l) by Pg, a, (3, and 7 , respectively, 
then we obtain dH). ■ 

Similarly, following the preceding derivation it can be shown that the eigenvalue distribution 
of interference satisfies (l22l) in the Stieltjes transform with ak, oi, 02 and 0:3 are replaced by Pj, 
K{L — 1)/M, P/M, and K{L — 1)/N, respectively. When the i.i.d. assumption on the channel 
SHi holds we have 

Corollary 2. Let M, K, and N are large with a, and y fixed yet not very small. Under rich scat¬ 
tering propagation (i.e., P 00 ), the eigenvalue distribution of F = Ps{SHiXi)\SHiXi)/MN 
converges to a non-random distribution with Stieltjes transform Gp satisfying 


7(1 + sGp) — PgGpi^ + sGp — 7 )(o + oisGp — 7 ) = 0. (23) 

Remark 2. Equation (1231) exactly reproduces the result in /l^ Eq.(75)], where it was derived 
under uncorrelated channel model and the convention that Pg = N/K. 

Corollary 3. Assume the channel model ([5]). Purther, assume K, P and N are large with 
K/P —)■ o' and 7 fixed but not very small. Then adding more antennas at the BS, the eigen¬ 
value distribution of F = Pg{SHiXiy{SHiXi)/MN converges to a non-random distribution 
resembling the eigenvalue distribution under uncorrelated channel with P receive antennas. In 
Stieltjes domain, the AED of F satisfies (1231) with a replaced by a'. 

Remark 3. From Corollaries \2} and |2] we observe that the gap between signal subspace and 
interference subspace under a physical channel with P AoAs and unlimited number of receive 
antennas M is equivalent to the gap under uncorrelated channel with M = P antennas at the 
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BS. Since blind subspace method relies entirely on the distance between the two subspaces, it 
turns out that its performance depends mainly on P rather than M. Nevertheless, the array 
gain due to the increase of M can still be partly retained, especially when signal subspace is 
detached from interference subspace. After that, the performance will be limited by the degrees 
of freedom of the channel, and hence it saturates. 

3) The support of the distribution 

In this part we will characterize the support of the distribution, which is important for knowing 
under what values of the system’s parameters the distribution yields two distinct eigenvalues 
clusters, and hence subspace separability. It is known that the Stieltjes transform is increasing 
on intervals on the real line outside the support of its distribution function lIT^ . Further, its inverse 
function s(x)0 is also increasing in these intervals only. Therefore, the endpoints of the support 
can be determined by finding the local extrema of the inverse function s{x). Unfortunately, these 
local extrema are notoriously difficult to express in a closed-form solution. So this paper is only 
devoted to deriving approximate formula without explicitly obtaining the extreme points. To that 
end, we plot s{x) for real x to find the regions where s{x) in increasing, and hence the support 
of the distribution, denoted S can be defined as the finite union of regions on the real line where 
s(a:) is not increasing, that is 



(24) 


x-^lx2 

Xl<X2 


To find the support of the distribution, it seems more adequate to use the h ^ h matrix G 2 
(i.e., a matrix of much smaller dimensions) rather than the n x n matrix Gi, since we are 
only interested in the nonzero eigenvalues. By the aid of (1211) one can verify that the Stieltjes 
transform Gg 2 corresponding to the S'-transform (fT9l) satisfies 


sGg2 + o,kGG2i—sGG2 + — — 1)(—sGgj + — — 1)x 
m m n n 

{-sGg2 + - - 1 ) + 1 = 0 . 

p n 


n 


(25) 


''The inverse function s(a;) is obtained as a solution for variable s in the fixed-point equation satisfying the Stieltjes transform, 
where a; is a real dummy variable. 
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By expanding (l25l) and after simple mathematieal manipulations we have 

Qkllx^s^ — aJlx^ (—34 + M + N + P) 

+ X {MNP + akhx {3ll + MN + {M + N) P - 2lk {M + N + P))) s (26) 

+ MNP - akx (-4 + M) {Ik - N) (4 - P) = 0 

where we have replaeed Gc 2 by the dummy variable x and m, n, and p by M, N, and 

P, respeetively. Note that (|2^ is a cubic equation in s, and its three roots, denoted s{x) = 

{si(a;), § 2 ( 3 ^), 53 ( 0 ;)} must be computed, and hence the support is defined by (l2^ . We remark 
that (|2^ can be also used to find the support of eigenvalue distribution of interference by noting 
that ttk = Pi and Ik = K{L — 1). 

4) Double-sided spectral analysis 

So far we have treated the signal subspace and the interference subspace separately. However, 
as far as systems with finite dimensions are concerned in practice, the two sets of eigenvectors 
corresponding to the signal and the interference are somehow interconnected. Thus the eigen¬ 
value distribution of the signal and the interference, should be studied jointly. For the sake of 
mathematical tractability, we again assume high SNR regime, i.e., we set W) = 0 in dS]). Let 


Gi = {SH)D^/‘^XX^D^/^{SH)yMN. (27) 

Instead, it is more convenient to study the eigenvalue distribution of 

G 2 = D^/'^XX^D^/‘^{SH)\SH)/MN (28) 

without accessing to the information about the eigenvalues of Gi. Further, let G2 = G21G22, 
where G 21 and G 22 are defined as follows: 

G 21 = D^^^XX^D^/^N, (29) 

G22 = {SHy{SH)/M ( 30 ) 


Remark 4 . Equation ([ 28 ]) is indeed intuitive in the sense that it gives an insight into the behavior 
of the gap between the signal subspace and the interference subspace as the ratio Pg / Pi varies. 
To see this, assume a rich scattering environment, i.e., P —)■ cxo. From central limit theorem, as 
N, M ^ 00 with K,L fixed, XX^N, and {SHy{SH)/M reduce to identity matrices. Thus 
the final spectrum of (l 28 l) is dictated by the spectrum of D, that is, a probability mass function 
(pmf) with two masses at Pg and Pi. In this case, subspace separability is possible whenever 
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Ps > Pj. However, when the dimensions of random matrices in (1281) grow large, but with fixed 
ratios the eigenvalue distribution converges to a non-random distribution rather than the pmf. 
This leads to eigenvalues spread around the center eigenvalue^ Therefore, a large difference 
between those deterministic parameters of the system is important which help to shift the two 
clusters of eigenvalues away from each other. 

From II 20 II . the asymptotic Stieltjes transform of X^DX /N, denoted Gxt£»x/Ai('S)’ 

N ^ 00 with KL/N ^ y is the unique solution of 

Xf{X)dX 




-1 


(31) 


1 + xi DX/N^^\ 

where /(A) is the probability density function (pdf) of eigenvalue of X^DX/N. It would seem 
difficult to evaluate (l3TI) . Since the entries of X are assumed i.i.d. CA/”(0,1), this class of matrices 
is unitarily invariant and thus asymptotically free w.r.t. any Hermitian matrix (see Appendix). 
Then it follows from |I21 Eq.(43)] that 

1 




-Sn{z) 


(32) 


yz + l 

which is derived within the framework of FPT. Since the eigenvalue pdf ff){x) of D converges 
to a pmf of two distinct eigenvalues, namely, Pg and Pj, hence ff){x) can be written as 


foi^) = - Ps) + - Pi) 


(33) 


where 5(.) is the Dirac delta function. By the definition of Stieltjes transform 0, [1121 . [[T3l . 
one can verify that ff,{x) admits a Stieltjes transform G^(s) of the form 

LPg — Ls + Pj — Pg 


L{Pg - s){Pi - s) 


(34) 


and from lfT2l Def. 3.4], it is easy to check that the corresponding S'-transform fulfils the following 
quadratic equation 

LPiPgzS^^iz) - {Pg - Pi + LPi + LPiz + LPgz)Sp,{z) + {L + Lz) = 0 (35) 


with the two roots 




b{z) ± ^Jh^iz) - AL‘^PiPg{z + l)z 


2LPjPgZ 


(36) 


^The center eigenvalues are determined by the deterministic quantities in the channel model. For instance, YiY^/MN has 
two deterministic quantities, namely, Ps and Pj corresponding to the effective signal and interference powers, respectively. 
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where b{z) = Ps — Pi + LPi + LPjz + LPgZ. One may prove that the root with minus sign 
is the true solution to ([35]^. Let KL/M —)■ a and KL/P t] as KL, M, and P ^ oo. 

Plugging (l3^ into (l32l) . combining the result with (fTVl) . and assuming the factors (l29l) and (l30l) 
are asymptotically free (see Appendix), then it follows from [fT3l Thm. 2.32] the S'-transform 
corresponding to G 2 can be finally written as 

b{z) - ^yb^z) -4L^PiPs{z + l)z 


2LPiPsz{'yz + l){az + l){riz + 1)' 

By making use of (|2TI) in (l37l) we have the following result. 


(37) 


Proposition 4 . Let S, H, D, and X be defined as in (|7]). Let also KL, N, P, and M ^ 00, 
with KL/M —)■ a, KL/P —)■ p, and KL/N 7. Further, denote the Stieltjes transform 
corresponding to G2 (l 2 ^ . Then Gq^ is the unique solution of the following fixed-point equation: 


LPs{sGq^ + 1) — Ps + P/(l + L{sGq^ + 2GPs{o.sGq^ 

+ a- IfipsG^^ +p- l){psG^^ + 7 - 1 ))) 

+ {PKLsGq^ + 1)^ + P^{LsGq^ + L - ly 

- 2PiPfil -L + L\sG^^ + l)sG^^)f/^ = 0. (38) 


Using (1381) to find the support of the distribution is difficult, hence, we adopt an approximate 
solution. Rewriting (IMI) in terms of the parameters K, L, M, N, and P yields 
2K^L^PiPs 


MNP 


- 2K‘^L^PjPM 


1 1 1 N 2 

■MN MP NP’ 


+ {2KL‘^PiPsx{— + jY + p) + + Ps))v 

+ ipmv-i)+ir+p!iLv-ir 

+ 2PiPfiL-L\v-l)v-l))P^ 


LPj + Pi — Pg — 2LPiPsX = 0 


(39) 


where we have replaced Gq^ by the dummy variable x and v = sx 1. Note that (l39l) (from 
which the all roots s{x) must be found and all extreme values are identified) is a polynomial of 
degree 6 in variables s and x, hence, finding its zeros is indeed intractable. Instead, we adopt 
an approximate solution such that the above higher order polynomial is reduced to a low-order 


®Since the S'-transform of a distribution with a single mass Ps is l/Ps, when Pj = Ps, it follows that the root with minus 
sign yields S-transform of the form 1/Ps- 


April 23, 2015 


DRAFT 










18 


polynomial of degree 2 in variable s. To that end, notiee that the multiplieative coefficients of 
the terms and v scale like a'q'y, aj + arj + rj'y, and a + rj + respectively. It turns out 

that in order to truncate out the terms and and hence have a good approximation of the 
true solution, the following conditions should be satisfied: 


a + V + 1 ^ ^ Q + ^ + 7 ^ ^ 


(40) 


arj'y a-j + ar] + tj'J 

where ignoring the terms and don’t affect the final answer significantly. Therefore, after 
truncating the higher order terms, namely, and we have 

{2KL^PiPsx{— + — + —) + L {Pi + Ps))v 

+ {P^{l + L{v-l)f + P^,{Lv-lf 

+ 2P/P,(L - L‘^{v - l)v - 

- LPi + Pi-Ps- 2LPiPsX = 0 (41) 

and its two roots si(x) and S 2 {x) can be computed, and thus the support is defined by (l24l) . 


B. Physical channel model with distinct AoAs 

So far we have considered the case of identical AoAs. However, depending on the complexity 
of the propagation environment, we may have different signals received from different/shared 
directions, captured by the channel model ([3]). Indeed, the analytical analysis of this channel 
model is extremely intricate. For the sake of simplicity, there is no significant loss of generality 
incurred by using the channel model dH). Here we recall that the received signals of the i\h cell 
impinge upon the array of the intended BS from i.i.d. AoAs of cardinality Pj. We also assume 
that all sets of AoAs over all cells are mutually independent. Without loss of generality, we 
consider the worst case of power imbalance between the desired signal and the interfering signal 
as has been stated in Remarl^ Now we may rewrite (@1) as 

L 

Yi = s/PsSiHiXi + ^i Y, + Wi. (42) 

Studying the eigenvalue distribution of the signal model (l42l) is intractable even when Wi = 0. 
Therefore, a seemingly natural way is to do one-sided analysis of the distribution of the first and 
second terms in (l42l) . This, however, would help us get some insights into the behavior of the 
eigenvalue distribution under the channel model (l42l) as it will be shown shortly. For the sake 
of illustration, we let / = 1. 
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First, note that the support of the eigenvalue distribution of the first term in (@21) can be 
computed by (l26l) . On the other hand, the second term of (l42l) can be rewritten in a compact 
way: 

Yi = (43) 


where Sj = {S 2 , S 3 , ■ ■ ■, Sl) G is the composite steering matrix, where n = P 2 +P 3 -\ -h 

Pl- Hj G is a block diagonal matrix whose blocks are H 2 , H 3 , ■ ■ ■, Hl and Xj = 

comprises all interference signals. Note that the eigenvalues of Hi 
equals the combined eigenvalues of the submatrices on its diagonal. Moreover, the eigenvalue 
distribution of each Hermitian matrix HiH} converges to a non-random distribution as K, Pi —)■ 
cx) with a fixed ratio K/ Pi. This implies that the eigenvalue distribution of HjHI also converges 
to a non-random distribution. To proceed with the analysis, the following Lemma enables us to 
characterize the behavior of eigenvalue distribution asymptotically. 

Lemma 5. Assume K < Pi,i = 2, ■ ■ ■, L with fixed ratios K/Pi ^ fii. Further, let A* = Pi/n and 
fi{x) be the eigenvalue distribution of HiHl- Then, as K,Pi — 00, the eigenvalue distribution 
fix) of HiH\ is related to fiix),i = 2, 3, • • •, L through the following relationship: 

L 

= (44) 

i=2 

Proof: Since the eigenvalues of Hj are equal to the combined eigenvalues of H 2 , ■ ■ ■ ,Hl, 
we can write 

triHiHj) = triH2Hl) + • ■ ■ + triHLHl) 


and the kth moment of the eigenvalue distribution of HjH] is given by 


triHiH])^ = triH2Hl)^ + • • ■ + triHLH 


Uk 


rUk 


(45) 


Note that the number of eigenvalues of HjHj are n, while each submatrix HiHj has Pi 
eigenvalues (i.e., the assumption K < Pi implies K nonzero eigenvalues and Pi — K zero 
eigenvalues). Thus the normalized trace form of (l45l) reads 


trniHiH\)^ = ^P2trpfiH2Hl)^ + ■ ■ ■ + PLtrp^iHLHl)^). 


Uk 




n 


(46) 


Recall that the Stieltjes transform G£(s) of n x n Hermitian matrix E can be expanded in a 
Laurent series involving the moments of E as (see e.g., tra Thm. 3.3]) 




1 A tr^iE’^ 


s ^' s" 

k=0 


(47) 
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Combining (l47l) with (|4^ yields 


L 

i=2 

Then the eigenvalue distribution is reeonstrueted from by applying the inversion 

formula of Stieltjes transform lfT3]l . ifT^ to obtain (1^ . ■ 

Since each submatrix Hi consist of i.i.d. complex Gaussian entries, it follows from (HD) that 
the eigenvalue distribution of HjHj is a weighted sum of the Marcenko-Pasture laws [fT^ with 
different ratios /9j, i = 2, • • •, L. One particular case is when P 2 = ■ ■ ■ = Pl = P, ii follows that 
the eigenvalue distribution of HjHj reduces to a Marcenko-Pasture distribution with a single 
parameter (5 = K/P. 

Proposition 6. Consider the channel model (l4^ . Without loss of generality, assume Pi > P 2 > 

■ ■ ■ > Pl- Then the users in the Lth cell contribute most to the spreading of the eigenvalues 
of channel. Further, under any system parameters, the power ratio Pg /Pi should be adjusted 
according to the Lth cell, which is required to maintain specified link performance when subspace 
method is used at BS. 

Remark 5. We observe that applying the free multiplicative law to the channel model (l4^ 
is extremely difficult because of nonidentical dimensions of the matrices Hi,i = 1,2, •••,//. 
Therefore, we provide here a non-rigorous yet intuitive proof of Proposition^ Further, the fact 
that the users experiencing less channel scatters render the eigenvalue distribution getting wider 
is indeed very intuitive. Notice that the maximum gap between the two subspaces occurs when 
P ^ 00 (i.e., uncorrelated channel), i.e., K/P —)■ 0. In addition, from Lemma |5] the eigenvalue 
distribution of HiH] is a weighted sum of Marcenko-Pasture distributions with different ratios. 
Assume that all limiting ratios, K/Pi, are infinitely small except one ratio, say the last ratio, 
which is made large (i.e., the signals of the Lth cell are received at the desired BS from a very 
small number of AoAs). Thus the support of the resulting distribution is dominated by the Lth 
distribution. This leads to a wider spread of the eigenvalues of the interference. 

To get insight into how the case of distinct AoAs may shape the distribution, and thus change 
the support of interference eigenvalue distribution, we assume the same number of AoAs per each 
cell w.r.t. the desired BS. Again by using one-sided spectral analysis and by following similar 
steps in Sec JTTT-A] with the aid of Lemma |5l one can prove that the eigenvalue distribution of 
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the N X N matrix YjYj/MN reads in the Stieltjes domain as 

- NPi{l + sx){K -KL + N + Nsx){N + P - LP 
+ Nsx)x + M{KPPix{L - 1)2 - N{L - 1) x 

(P + Pix{K + P))(l + sx) + Pix{N + NsxY) = 0. (49) 

By computing the three roots of (l49l) . s(x) = {si(a:), S 2 (x), 53 ( 0 :)}, then the support is defined 
by (1241). 


IV. Numerical results 

In this section, we use finite-size scenarios to show some numerical results that verify the 
theoretical results obtained in the asymptotic limit. We also provide simulation results for the 
uncoded bit error rate (BER) and compare the performance of subspace-based channel estimation 
scheme lIHl under the physical channel model and the i.d. channel. The pilot-based channel 
estimation scheme is also shown. In our simulations, the BS antennas are sparsely spaced at 
twice the carrier wavelength unless explicitly stated otherwise. 

A. The support of AED 

In Figsl4^and|4b]we use (|2^ and (l4TI) . respectively, to show the approximated boundaries of 
the eigenvalue distribution. The number of receive antennas M = 400, P = 200 AoAs, L = 4 
cells, K = 5 users per cell, the coherence time N = 1000 symbol periodjl the desired signal 
power Ps = 0.1(—lOdB), whereas the interference power P/ = 0.025 (~ —16dB), and high 
SNR is assumed (i.e., W = 0). Further, the support of distribution of the desired signal and the 
interference can be read on the right vertical axes, while on the left vertical axes is the exact 
noise-free support obtained for an i.d. channel with M = 400. 

As expected, the eigenvalue clusters of the signal and the interference exhibit larger spread 
from both sides of the support, and hence the gap between the two clusters decreases significantly. 
In contrast, under the same setting, but i.d. channel, the eigenvalues show less spread and the 
gap between the two clusters is more pronounced. In practice, especially with noisy samples, 
it turns out that, realistic physical channel, relative to i.d. channel, incurs higher power ratio 


’a is taken to be a relatively large in order to highlight the effect of the number of AoAs and power ratio, rather than the 
effect of the number of measurements. 
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(a) approximate solution using d26b (b) approximate solution using (ED 

Fig. 4. The nonzero support of the AED of YiY^/M, M = 400, P = 200, N = 1000, K = 5, L = 4, = — lOdB, 

Pi = -16dB. 


Ps/Pi such that the spectrum splits into signal and interference eigenvalues in order to maintain 
a comparable BER. Otherwise, a high frequency reuse factor, for example, would be needed, 
which leads to poor use of spectrum resources 

The supported boundaries computed by (|2^ and (|4TI) are compared to the histogram of 
nonzero eigenvalues in FigjSl Obviously, the histogram is composed of two bulks of eigenvalues 
clustered around two eigenvalues, 25 and 100 (center eigenvalues). The solid and dashed vertical 
lines correspond to the endpoints approximated by (|4TI) and (l2^ . respectively. Note that the 
approximated boundaries from (l4TI) are almost in agreement with the boundaries of the two 
bulks obtained from the histogram of eigenvalues. 

Figj^ shows the superimposed probability densities of eigenvalue for the physical channel 
(solid line) and the i.d. channel (dashed line). The parameters N, K, L, Pg, and Pj are the same 
as those in FigjU In the physical channel, P = 100 identical AoAs for all users w.r.t. BSl, 
and the number of receive antennas M is set to 600. On the other hand, the number of receive 
antennas in the case of i.d. channel is set to the number of AoAs, that is M = 100. Note that the 
eigenvalue pdf of YiY^/M almost matches the eigenvalue pdf of the corresponding i.d. channel. 
Thus, this result implies that the gap between the two bulks of eigenvalues (corresponding to the 
signal subspace and the interference subspace) of the physical channel with P AoAs cannot be 
further improved beyond what can be achieved when the channel is i.d. with P receive antennas 
as M increases unboundedly (see CorjS]). Nevertheless, array gain can be still retained as M 
increases, but ultimately the performance will saturate as will be shown in the next subsection. 
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eigenvalue 


Fig. 5. Normalized histogram of nonzero eigenvalues of YiYi /M, M = 400, P — 200, N — 1000, Ff = 5, L = 4, 
Ps = — lOdB, Pi — — 16dB. The solid and dashed lines are the approximated boundaries obtained by (ED and l l26t . respectively. 



Fig. 6. Empirical pdfs of nonzero eigenvalues of YiYi /M for physical and i.d. channels, N = 1000, K = 5, L = A, 
Ps = -lOdB, Pi = -16dB. 


Finally, Figs ITal and ITb] show the normalized histogram of eigenvalue distributions for physieal 
ehannel model with distinet AoAs. The parameters N, K, L, Pg, and Pj are the same as those 
in Fig® In FigjT^ we independently generate four sets of equal number of AoAs (Pi = P 2 = 
P 3 = P 4 = 200) eorresponding to eells 1, 2, 3, and 4 w.r.t. BSl. On the other hand, in FiglTbl 
we fix the number of AoAs of eells 1, 2, and 3 to 200 (i.e., a large number of AoAs), whereas 
the number of AoAs eorresponding to eell 4 is set to 20 (i.e., very small number of AoAs). 
Compared with the histogram in FigjT^ it is elear that the fourth eell renders the interferenee 
bulk of eigenvalues (left bulk) more wider and henee the gap is less pronouneed. 
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eigenvalue eigenvalue 

(a) Pi = P 2 = Ps = P 4 = 200 (b) Pi = P 2 = Pa = 200, P 4 = 20 


Fig. 7. Normalized histogram of nonzero eigenvalues of YiY^/M with distinct AoAs, M = 400, N — 1000, K = 5, L = 4, 
= -lOdB, Pi = -16dB. 


B. Bit error rate 

To give an intuitive feel of the effeet of the physical channel on the performance of subspace 
method, in Fig [T] we show the uncoded BER versus the power ratio Pj/Ps- In our simulation 
we consider the uplink in a four-cell network with 5 users per cell and QPSK modulation 
scheme. The coherence time of the channel N = 400 symbol periods, all signals are received 
from identical AoAs with cardinality P = 200 and per-user SNR = —5dB. To show the effect 
of increasing the number of receive antennas M, we use different values: M = 200,400,600, 
where the antenna elements are critically-spaced. We consider full reuse of K orthogonal pilot 
sequences across all cells. In addition, we use linear zero-forcing (ZF) for channel estimation 
(after applying SVD) and matched filter (ME) for data detection. The performance comparison 
with the classical pilot-based channel estimator is also shown, in which ZF and ME receivers 
are used to estimate the channel and detect data, respectively. 

From Fig {H we notice that the physical channel incurs loss of performance as opposed to the 
i.d. channel. For instance, when each BS is equipped with 400 antennas and the number of AoAs 
is 200, to achieve the same BER of 10“^, the power ratio should be roughly doubled, compared 
with the case of i.d. channel. Also, subspace-based scheme outperforms the pilot-based scheme 
in all cases except at the very high interference level, which is unlikely in practical scenarios. 
Further, the performance of the subspace-based scheme improves gradually with increasing M. 
As expected, its performance under this physical channel becomes closer to its performance under 
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Fig. 8. BER versus the power ratio for distinct AoAs, Pi = P 2 = P 3 ~ 100, M = 400, K = 5, L — 4:, N = 400, per-user 
SNR = —5dB, ZF is used for channel estimation and MF for data detection. 


an i.d. channel with M = 200 receive antennas. However, for very small interference levels, 
it further improves. Actually, this verifies the fact that when M grows large, the eigenvalue 
distribution of this physical channel with P AoAs becomes identical to that of an i.d. channel 
with M = P receive antennas. It should be noted that it still benefits partly from array gain 
due to increasing M, and this array gain becomes more useful when the two subspaces start 
to detach (i.e., when two subspaces overlap, the estimation error of channel will get larger and 
thus have a dominant impact on the performance). After that point, the performance becomes 
dominated by the degrees of freedom of the channel only. 

In FigH] we show the BER when distinct AoAs are used per each K users in each cell w.r.t. 
BSl. To highlight the degradation of the performance due to the cell with the smallest number 
of AoAs, we fix the number of AoAs of cell 1, 2, and 3 to 100, while we vary the number of 
AoAs of the fourth cell, P 4 = 10,20,50,100. It is clear that when the subspace-based scheme 
is used, the performance improves as P 4 increases, whereas it is almost the same for all values 
of P 4 when the pilot-based scheme is used. In subspace mehod, the gradual improvement of 
the performance can be interpreted as a consequence of gradual compression of the eigenvalues 
cluster of interference when P 4 increases (see Fig 17^ and FiglTbl). However, when P 4 becomes 
sufficiently large, we observe a slight improvement of the performance (i.e., a saturation effect of 
the performance). This is because the performance becomes limited again by the other interfering 
cells associated with smaller number of AoAs. Further, it is expected that when P 4 decreases, 
the degradation in performance becomes even worse when the interfering power from the fourth 
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cell is comparable to the power of the desired users. 



Fig. 9. BER when N ~ KL under an i.d. channel, M — 400, K — 15, per-user SNR = OdB, L = 4, ZF is used for channel 
estimation and MF for data detection. 

Finally, in FigUlwe simulate the BER when the number of users of the network is comparable 
to the coherence time of the channel. In these scenarios the all-orthogonal pilot-based scheme 
doesn’t work. Instead, we compare its performance with the classical pilot-based scheme. We 
assume an i.d. channel, Ff = 15, L = 4 and the coherence time N varies around KL = 60, i.e., 
N = 30, 60,120. It can be observed from Figl^that in all scenarios, especially the critical cases 
when iV = 30 and N = 60, the subspace-based scheme outperforms the pilot-based scheme. 
Its performance also improves with increasing N while the performance of pilot-based scheme 
is almost the same for all values of N. This means that subspace-based scheme exhibits better 
flexibility than the classical pilot-based scheme and all-orthogonal pilot-based scheme. 

V. Conclusion 

Pilot contamination is usually considered as a performance bottleneck in massive MIMO 
systems, if the conventional linear channel estimation method is used. As a departure from 
the linear estimation framework, the subspace-based scheme proposed in (SI exhibits better 
flexibility under various scenarios. In addition, it offers the potential of completely eliminating 
the pilot contamination effect. However, previous works are based on the common assumption 
of independent channels, which would be violated by real-world channels where the number of 
scatterers and AoAs might be both limited. 

For the more sensible physical channel model, in this paper, we derived approximate analytic 
expressions in Stieltjes domain for the corresponding eigenvalue distributions in both scenarios 
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of identical and distinct AoAs among the cells. The results demonstrate that the physical channel 
will exhibit a larger spectral spread, thus resulting in a smaller gap between the two eigenvalue 
clusters, which correspond to signal and interference subspaces. Therefore, to obtain the same 
performance as of the i.d. channels, we must guarantee a larger power ratio between the intended 
users and the interfering users. In particular, for the scenarios of distinct number of AoAs, the 
required power difference is determined by the cell with the smallest number of AoAs. Moreover, 
it is shown that adding more antennas does not significantly affect the above spectral gap, even 
though, the array gain could be reaped until saturation. 

For future research, since the blind subspace-based channel estimation scheme relies primarily 
on the eigenvalues spread of the channel (non-Bayesian or data-driven method), rather than 
the statistical properties of the channel (Bayesian method), the exploration of the possibility 
of combining the two methods is of importance from the perspective of system performance, 
especially in the scenario of bounded angular spreads. 

Finally, the results in this paper demonstrate the significant impact of the physical channel with 
finite AoAs on the performance of massive MIMO system, especially for the blind subspace- 
based channel estimation methods. 


Appendix 

DISCUSSION ON THE FREENESS CONDITION 

In this appendix, we prove the asymptotic freeness of the products of matrices in (fT^ 
and (|2^ . Let {Ei} G and {E-} G be two families of bi-unitarily invarianj^ ran¬ 
dom matrices, whose AED’s, as n —?■ oo, converge almost surely to non-random distributions 
with compacted supports. Further, let {Zj} G and {Z^} G be two families of 

non-random diagonal matrices with almost sure convergence of their AED’s to non-random 
distributions with bounded eigenvalues as p —)■ oo. Then from lfT9l Thm. 4.3.11], the family 
{.eI}, {Zj}, {Zj}}, G I,j,j G J is asymptotically free. It is important to note that 
it is not necessary for E* and E- to be square matrices. This can be seen from the fact that Zj 
and Zj can have different dimensions, e.g., an arbitrary number of the last entries can be zeros. 
Equivalently, Ej and E- can be considered as non-square matrices whereas Zj and Zj are two 
square matrices, such that the combination of all matrices is defined. 


bi-unitarily invariant matrix is a matrix whose the joint distribution of its entries is invariant when both left- and right- 
multiplied by unitary matrices. 
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Now, let us take the three cases: Ei = E- and Zj = Zj = I; Ei = E- and Z- = I-,Ei = E<. 
Then for each case we define the respective products: 

EiEj, i e I (50) 

E\{U^ZjU)E,, ielJeJ (51) 

ZjEiEjz], i & I, j,j & J (52) 

where the factor (UZiU^) in (ISTl) follows from the fact that Ei is bi-unitarily invariant and 
hence can be replaced by UEiV\ where U and V are unitary matrices. Note that ZjU is 
a unitarily invariant Hermitian matrix with AED converges to a non-random distribution (i.e., 
the eigenvectors are distributed in a maximally random way implying that ZjV fulfills the 
condition of |[I3 Thm. 4.3.11]), and by [HH Thm. 4.3.11] any subset of (l50l) . (ISTI) and (l52]) 
forms a family that is asymptotically free. 

Now, from (O, since the entries of Bk and Ck are assumed Gaussian with zero mean and unit 
variance, then they form bi-unitarily invariant random matrices each with the AED converges to 
the well-known Marcenko-Pasture law [fT^ . Then the Hermitian matrices (AkBkY(A^Bk) = 
and CkCl correspond, respectively, to (ISTI) and (l50]) . Eurther, because the as¬ 
sumption of i.i.d. of A\Ak incurs insignificant penalty compared with the Vandermonde matrix, 
it follows from (l50]) and (ISH) that B\, A\Ak, Bk, and CkCl form a family that is free 
asymptotically. 

Einally, the product D^/‘^XX'^D^/'^(SHy(SH) in (|2^ can be treated in the same way. 
Note that D^^'^XX^D^^'^ and (SHy(SH) are special instances of (l52]) and (ISTI) . respectively. 
Therefore, the asymptotic freeness property holds true for the products in (|2^ . 
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